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We compute the quasinormal frequencies of rotating black holes using the continued fraction 
method first proposed by Leaver. The main difference with former works, is that our results are 
obtained by a new numerical technique which avoids the use of two dimensional root-finding routines. 
The technique is applied to evaluate the angular eigenvalues of Teukolsky's angular equation. This 
method allows us to calculate both the slowly and the rapidly damped quasinormal frequencies with 
excellent accuracy. 
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I. INTRODUCTION 



Gravitational waves emitted by perturbed black holes are well known to be dominated, at late times, by quasinormal 
modes. They are characterised by complex frequencies with positive imaginary parts and represent therefore damped 
oscillations. The quasinormal frequencies of black holes in asymptotically flat spacetime are of great importance from 
the astrophysical point of view. They give us information about the main parameters of a black hole such as its mass 
and its angular momentum. It is for this reason that they have been extensively studied for more than forty years. 

The pioneering work on this topic was done by Regge and Wheeler [l| , who first studied linear perturbations of 
static black holes. They derived a second order ordinary differential equation which describes scalar, electromagnetic 
■ and gravitational perturbations of Schwarzschild black holes. The problem of separating the wave equation for Kerr 
rotating black holes, was solved by Teukolsky Q, who managed to obtain a couple of independent partial differential 
equations for the radial and angular part of the perturbation. An extensive study of black hole perturbation can be 
found in 0. 

Quasinormal modes were first found by Vishveshwara and by Press p| through numerical computations of the 
<«J ' time evolution of the gravitational waves around the blach hole. The first method for calculating them numerically 
', was proposed by Chandrasekhar and Detweiler [|| , and it was based on the direct numerical integration of the Regge- 
Wheeler equation. There have also been analytical attempts to obtain these frequencies and all the methods have also 
been applied to the Kerr and Reissner-Nordstrom black holes. However, these methods cannot provide very accurate 
values for the frequencies and break down for rapidly damped modes. 

Leaver , has shown that quasinormal frequencies of both static and rotating black holes can be determined by a 
continued fraction method. There are several methods to compute the quasinormal frequencies for static black holes, 
but the continued fraction method is still the only one which can be generalised to the Kerr case, where beside the 
quasinormal frequencies we must evaluate the separation constants, that is the angular eigenvalues associated with 
the Teukolsky equations. An improvement of Leaver's results can be found in 

In rotating black hole perturbation it is, indeed, the simultaneous computation of the angular eigenvalues that 
causes most of the problems of inaccuracy of the related frequencies. An alternative approach to Leaver's is to use 
ordinary perturbation theory, as Sasaki did to solve the angular equation, but this method only works as far as 
we consider small values for the expansion parameter aw. In other words, for non zero values of a, this technique 
does not allow us to deal with high frequencies. 

In this paper we follow Leaver's continued fraction method, but after obtaining an analytical solution for the radial 
and angular equations, we do not solve the related two continued fraction equations simultaneously as he did. Our 
approach consists, first of all, in the computation of the angular eigenvalues by a Pade' approximation of the true 
values, and then in the use of this approximation into the radial solution. 

Where previous results existed, our results are in close agreement and therefore provide independent confirmation. 
We also give new results for higher frequency modes and tabulate the Pade' approximations for the angular modes, 
which can be used in a wide range of applications concerning rotating systems. 
In what follows, we use geometric units, that is c = G = 1. 
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II. ANGULAR AND RADIAL CONTINUED FRACTION EQUATIONS 

In the Boyer-Lindquist coordinates (i, r, 6, <£>), rescaling t and r such that c = G = 1, the Kerr metric is 

, / 2Mr\ 9 / AMarsin 2 9\ S , , 9/99 2Ma 2 rsin 2 6>\ , , . 
(is 2 =11 — Ui 2 + - — tfr 2 - Ed6» 2 - sin 2 6 r 2 + a 2 + d<j) 2 . (2.1) 



V B / A 

Here M is the mass of the black hole, a its angular momentu m p er unit mass (0 < a < M), B = r 2 + a 2 cos 2 9, and 
A = r 2 — 2Mr + a 2 . The event horizon is at r = r+. Thorne [lQ| has shown that a ~ 0.9984M is probably an upper 
limit to the rotation realisable by an accreting astrophysical black hole. 

Scalar, electromagnetic and gravitational perturbations of this metric are all described by a master perturbation 
equation which, assuming a t and <fi dependence given by e -^ t + lm <f> ^ was separated by Teukolsky 0. 

The separated differential equation for the angular part of the perturbation is 

(l-« 2 )$T-2«^+W(«)5, = 0, (2.2) 



where 



and that for the radial part is 



229^ > (m + su) 2 

i—vr 



(2.3) 



where 



with 



A^i + 2(s + l)(r - M)^ + V(r)R s = 0, (2.4) 
ar z ar 



Tr K 2 — 2isK(r - M) A . 2 2 

V= ^ + Aisur - A + 2au)m - a u) , (2.5) 



u = cos9, (2.6) 
K = (r 2 + a 2 )uj - am. (2.7) 

The field spin-weight parameter s takes the values 0, — 1, —2, respectively, for outgoing scalar, electromagnetic, 
and gravitational waves. A is the angular separation constant for (2.2), and it reduces to 1(1 + 1) — s(s + 1) at the 
Schwarzschild limit. 

For rotating black holes, we have to solve (2.2) numerically following Leaver 0- Boundary conditions for (2.2) are 
that S s is regular at the regular singular points u — ±1. The indices there are given by ±(m + s)/2 at u = 1 and 
±(m -s)/2atit= — 1. A solution to (2.2) may be written as 

00 

S s (u) = e auJU {l + u)\ m ~ s ^ 2 {l - u)\ m+s \' 2 fin(l + u) n . (2.8) 

71=0 

The series coefficients are related by a three term recurrence relation and the boundary condition at it = 1 is satisfied 
only by its minimal solution sequence. The recurrence relation is 

a »i + Pqciq = 0, (2.9) 

a„a n+1 + p n d n + 7„a„_i = (n > 1), (2-10) 

where k\ = \m — s\/2, k2 — \m + s\/2 and the coefficients of the recurrence relation are 

an = -2(n+l)(n + 2fci + l), (2.11) 
(3 n = n(n - 1) + 2n(ki + k 2 + 1 - 2aoj) - 2aoj(2ki + s + 1) 

+(jfei + k 2 ){k l + k 2 + 1) - (a 2 uj 2 + s(s + 1) + A), (2.12) 

7„ = 2auj(n + ki + k 2 + s). (2.13) 
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The solution will be minimal if the angular separation constant A is a root of the continued fraction equation 

= 8 Q'oTi "172 «27 3 ^374 (2 l&\ 

0i- h- h-h-'"' 

A solution to the radial equation (2.4) is found similarly to the solution for (2.2). This equation has two regular 
singular points r + and r_ as the roots of A. We also define a new rotation parameter b = (1 — a 2 /M 2 )? , so that b 
varies from 1 to as a varies from to M (Kerr limit). Let us set a = a/2M and Co = 2Mu>. The event horizon is at 
r = r + . The indices at r — r + are ia + and — s — icr+, where er + = (tor + — am)/b. Since the second index corresponds 
to in-going radiation, we can establish on R s the quasinormal boundary conditions 

R s ~{r-r+)- s - l °+, asr^r+, (2.15) 
R s ~ r -l-2s+*v« ;T ", asr^oo. (2.16) 

Thus, our solution can be expressed as 

oo / v n 

R s = e iur (r - r _)-i— f^+«°-+( r _ r + )- s - ia + ^ a n ( - - r+ j , (2.17) 

where the expansion coefficients are again defined by a three term recurrence relation 

a a! + /Vo = 0, (2.18) 



The recursion coefficients are 



a n a n+1 + (3 n a n + 7„a„_i =0 (n > 1). (2-19) 

a„ = n 2 + (c + l)n + c , 

P n = -2n 2 + ( Cl +2)n + c 3 , 

7« = n 2 + (c 2 -3)n + c 4 - c 2 + 2, (2.20) 

and the intermediate constants are given by 

2i fCj „ 

cn = 1 — s — iio am 

b \2 

ci = -4 + 2icj(2 + b) + y ^~ - am j . 
C2 = s + 3 — ozcij — — I — — am 



c 3 = u) (4 + 2b - a) - 2amCo - s - 1 + (2 + &)wl> -AH ( — - am 



ACo + 2i ( Co 



c 4 = s + 1 - 2t2T - (2s + 3)iu I - - amj , (2.21) 

The radial series solution converges and the r = oo boundary condition is satisfied if, for a given a, m, A, and s, the 
frequency w is a root of the continued fraction equation 

= /30 "^^^'"' (2 ' 22) 

or any of its inversions. 

It is useful to bear in mind the behaviour of the quasinormal frequencies of Kerr black holes under complex 
conjugation. Looking at equations (2.2) and (2.4) we can verify that, if p nm = —iu> n ,m an d A; m are a quasinormal 
frequency and the corresponding angular eigenvalue, then p n ,- m — Pn.m an d A/._ m = A* m . This satisfies the 
requirement of complex conjugate quasinormal modes, because the perturbative equation was separated by Teukolsky 
2] , writing the wave equation as 

*(i,r,M) = ^J e- iwt J2 e im *S s (u)R s (r)cLj, (2.23) 

i=|s| rn— — l 

and in this equation the sum is over both positive and negative values of m. 
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III. NEW TECHNIQUE FOR THE ANGULAR EQUATION 



The angular separation constant can be approximated for small aw by perturbation theory. This is, for example, 
the method used by Sasaki in 9]. The angular eigenvalues are, in that case, expanded in powers of aco, 



A = A + acu\ 1 + a z u 2 \ 2 + 0{{au>) 6 ) 



(3.1) 



Explicit values of the coefficients are 



A„ 
Ai 



1(1+1) -2= (l-l)(l 
1(1 + 1) + 4 



= — 2m 



1(1 + 1) 



2a+i)(c!+ 1 ) 2 +2Uc;- i ) 2 + --- 



2 (I + 4)(l - 3)(l 2 + I - 3m 2 ) 

3 1(1 + l)(2l + 3)(2l - 1) ' 



with 



J+i 

-Irn 



(l + 3)(l-l)(l + m + l)(l - m + 1) 



"I 1/2 



c 1 - 1 - - 

c lm — 



(2l+l)(2l + 3) 



(I + 2)(l - 2)(l + m)(l - m) 
(2Z + 1)(2Z-1) 



1/2 



This method works as far as the expansion parameter aw is sufficiently small, so that if removes the possibility of 
dealing with high frequencies. Since we are interested both in slowly and rapidly damped modes, we prefer not to 
have such a limitation. As already mentioned in the introduction, we follow Leaver's continued fraction method and 
we adopt a new strategy for the numerical solution of the two continued fraction equations (2.14) and (2.22). Instead 
of a two dimensional root finding routine, which gives simultaneously frequencies and separation constants in the 
complex plane, we first compute the angular eigenvalues and then use those values in the continued fraction coming 
out from the radial solution. 

Firstly, we write a computer program able to approximate the continued fraction (2.14) in the most efficient way 
and calculate the roots, which are the angular eigenvalues. We use an excellent method known as modified Lent's 
algorithm |llj |. which essentially converts the continued fraction to a quotient of two series 



fn 



(3.2) 



where A n and B n are determined from recurrence relations connected to the the coefficients of the continued fraction. 
The series are evaluated from left to right, and the algorithm is stopped when the value of the absolute value of the 
difference between two consecutive fractions is sufficiently small. Within the same program, the roots are searched 
by the use of the complex one dimensional Ncwton-Rapson method routine. 

Then, we fit these numerical data for A to a Pade' approximation for real atu given by 



A 



a 



- aiaui + a2(aui) 2 + a^au) 3 + a^aui) 4 
1 + b\auj + b2(auj) 2 + ^(aw) 3 



(3.3) 



Finally, the Pade' approximation and numerical data are compared for imaginary auj. 

In the Pade' formula (3.3), the numerator is one degree higher than the denominator. The reason for this choise is 
that, the Teukolsky's angular equation belongs to the family of differential equations for angular prolate spheroidal 
wave functions. The behaviour of the angular eigenvalues, for large ato, related to this kind of equation is well known 

m 

Fitting simulations have been done for different values of the angular parameter a, for several values of the spherical 
harmonic parameters I and m and they always work accurately when we compare our results with Leaver's. In Fig.l, 
we show how well our Pade' approximation works for m = and I = 2 both for the real and imaginary parts of A. 
The fitting routine we have used is gnufit, but the same results can be obtained using nonlinear curve fitting packages 
in Maple, available on the web. The fit parameters of the rational function representing A are shown in Table 1 for 
m = and in Table 2 and Table 3 for m = ±1. 



■5 




FIG. 1: The real and imaginary parts of the Pade' approximation of the angular eigenvalue A are shown on the left and right 
columns respectively. They are plotted as functions of lu. These figures represent results we get for I — 2, m = and with the 
angular parameter a = 0.1, 0.3, 0.49 from the top couple of plots to the bottom one. Note that, the Pade' fit of the numerical 
data for A has been done for real values of w. 



IV. RADIAL EQUATION AND NUMERICAL RESULTS 

As a consequence of the new technique, the angular equation with its angular eigenvalues is not an issue anymore. 
We are in the position of carrying on our numerical search for quasinormal frequencies and we switch all our attention 
at the radial equation. 

As already done in the angular case, we must efficiently represent the radial continued fraction (2.22). It is 



6 



a 


a 


cii 


ai 


a3 


«4 


bi 


b 2 


bi 


0.0 


4.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.1 


4.0 


-0.2979 


-0.3994 


0.0420 


-0.0294 


-0.0744 0.0310 


0.0009 


0.2 


4.0 


-0.3127 


-0.3939 


0.0434 


-0.03 


-0.0781 


0.0323 


0.0008 


0.3 


4.0 


-0.5886 


-0.2857 0.0679 


-0.0396 


-0.1470 


0.0588 


-0.0009 


0.4 


4.0 


-0.5284 


-0.2889 


0.0619 
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0.0627 


-0.0404 
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0.0586 


-0.00086 



TABLE I: This table shows, for m = 0, the parameters we obtain by fitting the Pade' approximation for A, to the numerical 
data coming out from the computer program for the calculation of the angular eigenvalues of Teukolsky equation. 
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TABLE II: This table shows, for m = —1, the parameters we obtain by fitting the Pade' approximation for A, to the numerical 
data coming out from the computer program for the calculation of the angular eigenvalues of Teukolsky equation. 
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TABLE III: This table shows, for m — +1, the parameters we obtain by fitting the Pade' approximation for A, to the numerical 
data coming out from the computer program for the calculation of the angular eigenvalues of Teukolsky equation. 



convenient in this case if we consider the 71th inversion of the radial continued fraction, that is 



= 



Pn 



Otn-lln a n -27n-l 



"071 



Pn- 



Pn- 



Ctnln+1 a n +l"/n+2 



Pn 



Pn 



+2" 



(4.1) 



because the nth quasinormal frequency is usually found to be numerically the most stable root of the nth inversion 
0. We, then, use again Lentz's method both for the finite and the infinite fractions in (4.1). Finally, we write another 
simple computer program similar to the one used for the angular equation, able to approximate the sum (4.1) and at 
the same time to compute its roots. In this second program we make use of the Pade' rational function to represent 
the angular eigenvalues which appear in (2.22). In this case the roots are obviously the Kerr quasinormal frequencies 
we are looking for. 

Using this new approach we actually recover former results up to the ninth mode 0,0 and in some cases we are able 
to go even beyond them. In Fig. 2, we have some plots showing patterns of quasinormal frequencies for the angular 
parameter a ranging from a — 0.1 to a = 0.49, which is very close to the Kerr limit. The top left box represents the 
frequencies for a = 0.0, I = 2, and m = and it is the classical result for Schwrzschild black holes. It can be found, for 
instance, in ^ n Fig- 3, we show details of the modes on a different scale. Parametrisation is done by the the 

value of the rotation parameter a. It looks clear that the ninth mode is initially on the imaginary axis as expected. 

In Fig. 4, we plot a detail of the curve representing the fifth mode for m = ±1 and we illustrate the complex 
conjugate symmetry of the Kerr quasinormal frequencies mentioned at the end of section 2. This figure, shows again, 
that our technique reproduces Leaver's results accurately. The peak we get for m = — 1 and angular parameter a 
ranging between zero and the Kerr limit, confirms similar behaviour found in Leaver 0. 

In Table 3, we give the numerical values of the Kerr quasinormal frequencies corresponding to I = 2 and m = 0, 



OJR LOR 

FIG. 2: First quasinormal frequencies for I = 2 and m = scaled by 2M, represented in the complex plane. In each box we 
have a different value for the angular momentum per unit mass a, which ranges between 0.0 and M (very close to the Kerr 
limit). 



obtained with the new technique. Note that, the values for the fundamental modes can be compared with Leaver's 
results and they match accurately, see Table 4. 
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FIG. 3: The first, third, fifth, seventh and ninth quasinormal modes for I = 2 and m — 0, scaled by 2M, are shown. The single 
modes are paramtrised by the value of the rotation parameter a. The left upper end of each mode (squares), represents the 
frequency at the Schwarzschild limit a = 0. The right end in the first and third modes and the lower left end in the others, 
represent the frequency at a = 0.49. The intermediate dots, correspond to a = 0.1,0.2,0.3,0.4. As clearly shown, the ninth 
mode is on the imaginary axis for a = 0.0, after that it moves away as expected. 



(0.747343, -0.177925) 
(0.750252, -0.177401) 
(0.75936, -0.175653) 
(0.776104, -0.17199) 
(0.803834, -0.164314) 
(0.844509, -0.147065) 



(0.693422, -0.547829) (0.602107, -0.956555) 

(0.697296, -0.546029) (0.607671, -0.952749) 

(0.709343, -0.540027) (0.624762, -0.940155) 

(0.731061, -0.527509) (0.654649, -0.914241) 
(0.765136, -0.501476) (0.6975, -0.86143) 
(0.802873, -0.446104) (0.7198, -0.7616) 



(0.503011, -1.4103) (0.415028, -1.89369) 

(0.510396, -1.40365) (0.424195, -1.88378) 

(0.532732, -1.38181) (0.451287, -1.85134) 

(0.57023, -1.33741) (0.494169, -1.78589) 

(0.616905, -1.24857) (0.535028, -1.65631) 

(0.60027, -1.11018) (0.462988, -1.50839) 



TABLE IV: In this table, we show the numerical values of the first five Kerr quasinormal frequencies for I 
a = 0.0, ..,0.49. 



2, m = and 



a 



Leaver's fund, mode Our fund, mode 



0.0 
0.1 
0.2 
0.3 
0.4 
0.49 



(0.747343, -0.177925) (0.747343, -0.177925) 

(0.750248, -0.177401) (0.750252, -0.177401) 

(0.759363, -0.175653) (0.75936, -0.175653) 

(0.766108, -0.171989) (0.776104, -0.17199) 

(0.803835, -0.164313) (0.803834, -0.164314) 

(0.844509, -0.147065) (0.844509, -0.147065) 



TABLE V: In this table, we compare our fundamental mode to Leaver's, for 
results are. 



2 and m — and we show how close the two 
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-1.76 



-0.6 



-0.4 



-0.2 



0.2 



0.4 



0.6 
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FIG. 4: A detail of the curve representing the fifth mode for m = ±1 is plotted. We also illustrate the complex conjugate 
nature of the Kerr quasinormal frequencies. The squares represent the complex conjugate frequencies for a = 0.0. 



V. CONCLUSION 



The quasinormal frequencies for the gravitational perturbation of a Kerr black hole are the most interesting from 
an astrophysical point of view and they have been matter of research for many years. 

Perturbations of rotating black holes can be reduced to the study of two continued fraction equations. Their 
common zeros are the angular and radial quasinormal frequencies. Searching for the solutions in a reliable way is a 
difficult task in numerical analysis. They could be found if we mapped out the full zero contours of both functions 
of the system. These zero contours would consist of an unknown number of disjoint closed curves. The roots are the 
intersections of pairs of zero contours • 

In this paper we have introduced a new technique to compute such frequencies in a way that avoids the use of 
four dimensional root finding routines. The new technique consists of a new numerical way to evaluate the angular 
eigenvalues of the Teukolsky angular equation, which is independent from the solution of the radial part of the problem. 

The results we have obtained work for any value of the spherical harmonic parameters I and m and for the angular 
parameter a ranging from the static to the Kerr limit. These results accurately agree with former work and results 
are given for higher modes than hitherto. It would be interesting to analyse the possibility to extend the numerical 
computation to higher modes or to look for an analytical formula to describe the asymptotic behaviour of the Kerr 
quasinormal frequencies in a way similar to what has been already done in the static case |13| . 

It is intended that the C+ + programme used to obtain the results in this paper, will be made available to the 
worldwide community as a useful tool in gravitational wave research. 
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